#=
测试置换
=#


include("../src/pimary.jl")
using ..pimary

using Plots
using LinearAlgebra


"""
绘制坐标
"""
function draw_sket!(plt, Nt, vol)
    for it in Base.OneTo(Nt)
        plot!(plt, [0.0, vol[1]], [it, it], line=:dot, label="", color=:black)
    end
    plot!(plt, [0.0, vol[1]], [Nt+1, Nt+1], line=:dot,
    xlims=(0.0, vol[1]), ylims=(0.0, Nt+2), label="", color=:black)
end


"""
绘制beads
"""
function draw_beads!(plt, ti, rpos, clr)
    for i in Base.OneTo(length(rpos))
        annotate!(plt, rpos[i], ti, ("$(i)", 10, :black, :center))
        #plot!(plt, Plots.partialcircle(0, 2π, 100))
    end
    scatter!(plt, rpos, repeat([ti], length(rpos)), markersize=15, label="", color=clr)
end


"""
绘制连线
"""
function draw_conns!(plt, ti, rpos, rprev, clr)
    for i in Base.OneTo(length(rpos))
        plot!(plt, [rpos[i], rprev[i]], [ti, ti-1], label="", color=clr)
    end
end



"""
测试置换求和==det
"""
function run1()
    melec = 1.0
    Δβ = 0.1
    hbar = 1.0
    ther = sqrt(2π*(hbar^2)*Δβ/melec)
    #
    Np = 3
    Nt = 1
    vol = [1.0, 1.0, 1.0]
    fbup = WrdFb(Float64, Np, Nt)
    uppos = rand(Np, 3).*vol
    for ip in Base.OneTo(Np); for it in Base.OneTo(Nt)
        fbup[ip][it, :] = uppos[ip, :]
    end; end
    println(fbup)
    r1 = beads(fbup, 1)
    r2 = copy(r1)
    #假设有一个小的位移
    for ip in Base.OneTo(Np)
        r1[ip, 1] += (rand()-0.5)*ther
        r1[ip, 2] += (rand()-0.5)*ther
        r1[ip, 3] += (rand()-0.5)*ther
    end
    plt = plot()
    draw_sket!(plt, Nt, vol)
    draw_beads!(plt, 1, r1[:, 1], :red)
    draw_beads!(plt, 2, r2[:, 1], :red)
    draw_conns!(plt, 2, r2[:, 1], r1[:, 1], :red)
    savefig("test.png")
    #三个粒子的交换一共6种情况
    #（1，2，3），（2，3，1），（3，1，2），（1，3，2），（3，2，1），（2，1，3）
    Z = 0
    #123
    #<r1r2r3| exp^{-π(rpp-rpos)^2/λe^2} |r'1r'2r'3>
    r2c = r2[[1,2,3], :]
    Z += 1.0 * free_kinetic(r1, r2c, vol; λe=ther)
    r2c = r2[[2,3,1], :]
    Z += 1.0 * free_kinetic(r1, r2c, vol; λe=ther)
    r2c = r2[[3,1,2], :]
    Z += 1.0 * free_kinetic(r1, r2c, vol; λe=ther)
    #
    r2c = r2[[1,3,2], :]
    Z -= 1.0 * free_kinetic(r1, r2c, vol; λe=ther)
    r2c = r2[[3,2,1], :]
    Z -= 1.0 * free_kinetic(r1, r2c, vol; λe=ther)
    r2c = r2[[2,1,3], :]
    Z -= 1.0 * free_kinetic(r1, r2c, vol; λe=ther)
    println(Z)
    #
    efk = exch_free_kinetic(r1, r2, vol; λe=ther)
    println(det(efk))
end


"""
测试对lndet进行求导的功能
"""
function run2()
    melec = 1.0
    Δβ = 0.1
    hbar = 1.0
    ther = sqrt(2π*(hbar^2)*Δβ/melec)
    #
    Np = 3
    Nt = 1
    vol = [1.0, 1.0, 1.0]
    fbup = WrdFb(Float64, Np, Nt)
    uppos = rand(Np, 3)
    for id in Base.OneTo(3)
        uppos[:, id] = uppos[:, id]*vol[id]
    end
    for ip in Base.OneTo(Np); for it in Base.OneTo(Nt)
        fbup[ip][it, :] = uppos[ip, :]
    end; end
    println(fbup)
    r1 = beads(fbup, 1)
    r2 = copy(r1)
    #假设有一个小的位移
    offs = rand(Np, 3)
    for ip in Base.OneTo(Np)
        r1[ip, 1] += (offs[ip, 1]-0.5)*ther
        r1[ip, 2] += (offs[ip, 2]-0.5)*ther
        r1[ip, 3] += (offs[ip, 3]-0.5)*ther
    end
    #交换矩阵
    mat = exch_free_kinetic(r1, r2, vol; λe=ther)
    println(mat)
    therprev = ther
    #对交换矩阵的矩阵元的求导
    melec = 1.0
    Δβ = 0.101
    hbar = 1.0
    ther = sqrt(2π*(hbar^2)*Δβ/melec)
    fbup = WrdFb(Float64, Np, Nt)
    #uppos = vol*rand(Np, 3)
    for ip in Base.OneTo(Np); for it in Base.OneTo(Nt)
        fbup[ip][it, :] = uppos[ip, :]
    end; end
    println(fbup)
    r1 = beads(fbup, 1)
    r2 = copy(r1)
    #假设有一个小的位移
    #offs = rand(Np, 3)
    for ip in Base.OneTo(Np)
        r1[ip, 1] += (offs[ip, 1]-0.5)*ther
        r1[ip, 2] += (offs[ip, 2]-0.5)*ther
        r1[ip, 3] += (offs[ip, 3]-0.5)*ther
    end
    mat2 = exch_free_kinetic(r1, r2, vol; λe=ther)
    println(mat2)
    #
    dmat = exch_derv_beta(r1, r2, vol; λe=therprev)
    println(dmat/0.1)
    println((mat2-mat)/0.001)
    #
    println(log(det(mat2)) - log(det(mat)))
    invmatT = inv(mat)
    dev = tr(dmat*invmatT)/0.1
    println(dev)
end

run1()
run2()

